function y = hdr_SampleCoords(hfile, c, vol, imeth, ctype)
% HDR::SampleCoords  - sample VoxelData at BV based coords
%
% FORMAT:       y = hdr.SampleCoords(c [, vol [, imeth [, ctype]])
%
% Input fields:
%
%       c           3xC or Cx3 coordinates
%       vol         volume number (default 1)
%       imeth       method 'linear', {'cubic'}, 'lanczos3', 'nearest'
%       ctype       coordinate type {'tal'}, 'bvsys', 'bvint'
%
% Output fields:
%
%       y           double precision data sampled at given coordinates

% Version:  v0.7f
% Build:    9030320
% Date:     Mar-03 2009, 8:52 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'hdr') || ...
   (~isa(c, 'double') && ...
    ~isa(c, 'single')) || ...
    ndims(c) > 2 || ...
   ~any(size(c) == 3) || ...
    any(isinf(c(:)) | isnan(c(:)) | abs(c(:)) > 256)
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
cfr = hdr_CoordinateFrame(hfile);
if numel(cfr.Dimensions) < 3 || ...
    ndims(bc.VoxelData) < 3 || ...
    size(bc.VoxelData, 1) ~= cfr.Dimensions(1) || ...
    size(bc.VoxelData, 2) ~= cfr.Dimensions(2) || ...
    size(bc.VoxelData, 3) ~= cfr.Dimensions(3)
    error( ...
        'BVQXfile:InvalidObject', ...
        'VoxelData and Dimensions mismatch.' ...
    );
end
if nargin < 3 || ...
    numel(vol) ~= 1 || ...
   ~isa(vol, 'double') || ...
    isinf(vol) || ...
    isnan(vol) || ...
    vol < 1 || ...
    (vol > 1 && ...
     (ndims(bc.VoxelData) < 4 || ...
      vol > size(bc.VoxelData, 4)))
    vol = 1;
else
    vol = fix(vol);
end
if nargin < 4 || ...
   ~ischar(imeth) || ...
   ~any(strcmpi(imeth(:)', {'cubic', 'lanczos3', 'linear', 'nearest'}))
    imeth = 'cubic';
else
    imeth = lower(imeth(:)');
end
if nargin < 5 || ...
   ~ischar(ctype) || ...
   ~any(strcmpi(ctype(:)', {'bvint', 'bvsys', 'tal'}))
    ctype = 'tal';
else
    ctype = lower(ctype(:)');
end
if size(c, 1) == 3 && ...
    size(c, 2) ~= 3
    c = c';
end
try
    tmat = inv(double(cfr.Trf));
catch
    error( ...
        'BVQXfile:InvalidArgument', ...
        'Coordinate transformation matrix not invertible.' ...
    );
end

% get coordinate system right
if strcmpi(ctype, 'bvint')
    c = 128 - c([3, 1, 2], :);
elseif strcmpi(ctype, 'bvsys')
    c = 128 - c;
end

% get sampling coordinates
c(:, 4) = 1;
c = double(tmat * c')';

% get data
if ndims(bc.VoxelData) > 3
    vdt = bc.VoxelData(:, :, :, vol);
else
    vdt = bc.VoxelData(:, :, :);
end

% scaling
if any([2, 4, 8, 130, 132, 136] == bc.ImgDim.DataType) && ...
   (bc.ImgDim.ScalingIntercept ~= 0 || ...
    bc.ImgDim.ScalingSlope ~= 1)
    vdt = bc.ImgDim.ScalingIntercept + bc.ImgDim.ScalingSlope .* double(vdt);
end

% sampling
offsamp = min(vdt(:));
offsamp = offsamp - offsamp;
y = flexinterpn_method(vdt, c(:, 1:3), offsamp, imeth);
